
##############################
### Income inequality data ###
##############################

rm(list = ls())
gc()

library(data.table)
library(ggplot2)
library(binsreg)
library(tidyverse)

### Set WD
setwd("~/Dropbox (Princeton)/Research/*Spatial_Inequality/paper_SMD_inequality/JoP_final/replication")  

### Load data
load("data/inequality.RData")

### Plots
pdf("out/fig3a.pdf", width=4, height=3.5)
(p <- (ggplot(data=inequality[ pcon_name == "uk" ], aes(y=diff_mean_med_tot, x=year))
       + geom_line(linewidth=0.9, show.legend = F, colour="black")
       + stat_summary(data=inequality[ is_region == 0 ], aes(y=diff_mean_med_tot, x=year), linetype="solid", geom="line", size=0.3, fun="median", colour="blue", show.legend = F)
       + stat_summary(data=inequality[ is_region == 0 ], aes(y=diff_mean_med_tot, x=year), linetype="dashed", geom="line", size=0.9, fun="median", colour="blue", show.legend = F)
       + annotate(geom = "text", x = 2017, y = 10500, label = "UK average", size=3.25, colour="black")
       + annotate(geom = "text", x = 2017.3, y = 7100, label = "Median constituency", colour="blue", size=3.25)
       + annotate(geom="text", x=2011.6, y=10700, label="more unequal", color="gray50", size=2.8)
       + annotate("segment", x = 2011, xend = 2011, y = 9000, yend = 10500, size=0.4, colour="gray50", arrow=arrow(length=unit(0.30,"cm"), type = "closed", angle = 20))
       + theme_bw()
       + ylab("Mean-median inequality (£)") 
       + xlab("Year")
       + theme(legend.position="bottom", panel.border = element_rect(color="gray"), panel.grid = element_blank()) 
       + guides(colour = guide_legend(nrow=2, byrow=T))
       + scale_x_continuous(breaks=seq(2011,2019,2))
))
dev.off()



### Binscatter
bs_plot <- binsreg::binsreg(inequality[ is_region == 0 ]$diff_mean_med_tot, inequality[ is_region == 0 ]$pop_density)

bs_out <- as.data.table(bs_plot$data.plot$`Group Full Sample`$data.dots)

setnames(bs_out, old=c("x","fit"), new=c("pop_density","diff_mean_med_tot"))

pdf("out/fig3b.pdf", width=4, height=3.5)
(plot <- (ggplot(data=inequality[ is_region == 0 ], aes(y=diff_mean_med_tot, x=pop_density))
          + geom_point(data=bs_out, alpha=0.5)
          + stat_smooth(method="loess")
          + theme_bw()
          + annotate(geom="text", x=1500, y=38000, label="more unequal", color="gray50", size=2.8)
          + annotate("segment", x = 0, xend = 0, y = 30000, yend = 37000, size=0.4, colour="gray50", arrow=arrow(length=unit(0.30,"cm"), type = "closed", angle = 20))
          + ylab("Constituency mean-median inequality (£)")
          + xlab("Population density")
          + coord_cartesian(ylim=c(5000,38000))
          + theme(legend.position="bottom", panel.border = element_rect(color="gray"), panel.grid = element_blank()) 
          + guides(colour = guide_legend(ncol=1, byrow=T))
))
dev.off()





### Summary statistics for inequality  
income_comb_sum <- inequality[ is_region == 0  ]

d1.1 <- income_comb_sum[ year == 2019 ,.(pcon,pcon_name,year,diff_mean_med_tot) ] %>% 
  arrange(desc(diff_mean_med_tot)) %>% 
  slice(1:10)

d1.2 <- income_comb_sum[ year == 2019 ,.(pcon,pcon_name,year,diff_mean_med_tot)] %>% 
  arrange(-desc(diff_mean_med_tot)) %>% 
  slice(1:10)

d2.1 <- income_comb_sum[ year == 2011 ,.(pcon,pcon_name,year,diff_mean_med_tot)] %>% 
  arrange(desc(diff_mean_med_tot)) %>% 
  slice(1:10)

d2.2 <- income_comb_sum[ year == 2011 ,.(pcon,pcon_name,year,diff_mean_med_tot)] %>% 
  arrange(-desc(diff_mean_med_tot)) %>% 
  slice(1:10)

d1x <- cbind(d1.1[ ,.(pcon_name,diff_mean_med_tot)],
             d1.2[ ,.(pcon_name,diff_mean_med_tot)])

names(d1x) <- c("Constituency_top","Inequality_top","Constituency_bot","Inequality_bot")

### SI Table B1
print(xtable::xtable(d1x, include.colnames=T, caption="Summary statistics", digits=0), include.rownames=F)


### SI Table B.2
d3.1 <- income_comb_sum[ year == 2019 ,.(pcon,pcon_name,year,diff_mean_med_tot,pop_density) ] %>% 
  arrange(desc(pop_density)) %>% 
  slice(1:10)

print(xtable::xtable(d3.1[ ,.(pcon_name,round(pop_density,0),diff_mean_med_tot)], include.colnames=T, caption="Summary statistics", digits=0), include.rownames=F)

### SI Table B.3
d3.2 <- income_comb_sum[ year == 2019 ,.(pcon,pcon_name,year,diff_mean_med_tot,pop_density)] %>% 
  arrange(-desc(pop_density)) %>% 
  slice(1:10)

print(xtable::xtable(d3.2[ ,.(pcon_name,round(pop_density,0),diff_mean_med_tot)], include.colnames=T, caption="Summary statistics", digits=0), include.rownames=F)



### Plots
ineq_v1 <- inequality[ is_region == 0 ]
setnames(ineq_v1, "pcon", "pcon_code")
setnames(ineq_v1, "diff_mean_med_tot", "diff_mean_med_tot_inc")


### UK wide inequality
ineq_UK <- inequality[ pcon_name == "uk" ,.(year, diff_mean_med_tot) ]
setnames(ineq_UK, "diff_mean_med_tot", "diff_mean_med_tot_inc_UK")

ineq_comb <- merge(ineq_v1, ineq_UK, by="year")

inc_2011 <- ineq_comb[ year == 2011,.(pcon_code,pcon_name,diff_mean_med_tot_inc)]
inc_2019 <- ineq_comb[ year == 2019,.(pcon_code,pcon_name,diff_mean_med_tot_inc,pop_density)]

setnames(inc_2011, "diff_mean_med_tot_inc", "diff_mean_med_tot_inc_2011")
setnames(inc_2019, "diff_mean_med_tot_inc", "diff_mean_med_tot_inc_2019")

inc_2011_19 <- merge(inc_2011, inc_2019, by=c("pcon_code","pcon_name"))

inc_2011_19[ , pct_chg := (diff_mean_med_tot_inc_2019-diff_mean_med_tot_inc_2011) / diff_mean_med_tot_inc_2011 *100 ]
inc_2011_19[ , pct_chg := as.numeric(pct_chg)]

pdf("out/si_fig_b1a.pdf", width=4, height=3.5) # uk_ineq_pcon_chg1.pdf
(plot <- (ggplot(data=inc_2011_19, aes(y=log(diff_mean_med_tot_inc_2019), x=log(diff_mean_med_tot_inc_2011)))
          + geom_abline(slope=1, colour="gray25")
          + geom_point(shape=21, fill="blue", size=1.5)
          + theme_bw()
          + coord_cartesian(ylim=c(7,12), xlim=c(7,12))
          + ylab("Inequality 2019, log")
          + xlab("Inequality 2011, log")
          + theme(legend.position="right", panel.border = element_rect(color="gray"), panel.grid = element_blank())
))
dev.off()


pdf("out/si_fig_b1b.pdf", width=4, height=3.5) # uk_ineq_pcon_chg3.pdf
(plot <- (ggplot(data=inc_2011_19, aes(x=pct_chg))
          + geom_density(fill="blue")
          + geom_vline(xintercept = 0, linetype="dashed", colour="gray50")
          + theme_bw()
          + ylab("Density")
          + xlab("Percent change, 2011-2019 (%)")
          + theme(legend.position="right", panel.border = element_rect(color="gray"), panel.grid = element_blank())
))
dev.off()






